
%% 1D Estimation.  Looking at Likelihood and Bayesian Posterior

% web('https://en.wikipedia.org/wiki/Bayes'_theorem')
% web('https://upload.wikimedia.org/math/1/4/a/14a536ae3b446d598d5d6091fee8298f.png')

A=1;
N=10;
X=0:N-1;
Sigma=3;

Expected=X*A;
Noisy=Expected+Sigma*randn(size(Expected));

figure
plot(X,Noisy,'o')
hold on
plot(X,Expected,'linewidth',2)
xlabel('x')
ylabel('y')
legend('Noisy Data','Expected Signal')

LogL=@(Theta,Data,Sigma,X)sum(-log(2*pi*Sigma^2) -(Data-(Theta*X)).^2./(2*Sigma^2));

ThetaTest=linspace(-1,3,200);
LogLResult=zeros(size(ThetaTest));
for nn=1:length(ThetaTest)
    LogLResult(nn)=LogL(ThetaTest(nn),Noisy,Sigma,X);
end

% From Bayes' rule
% P(Theta|Data)=P(Data|Theta)*P(Theta)/P(Data)

% This is a 'uniform' prior
PTheta = 1/length(ThetaTest)*ones(size(ThetaTest));

%This is the integral of the numerator in Bayes' rule. This acts as a
%normalization constant. 
DTheta=ThetaTest(2)-ThetaTest(1);
PData=sum(PTheta.*exp(LogLResult)*DTheta);

%This is the 'Posterior Distribution'
PTheta_Data=PTheta.*exp(LogLResult)/PData;


figure

%Show likelihood and prior
subplot(1,2,1)
[hAx]=plotyy(ThetaTest,exp(LogLResult),ThetaTest,PTheta);
xlabel('\Theta')
ylabel(hAx(1),'P(Data|\Theta)')
hold on
plot(ThetaTest,PTheta)
ylabel(hAx(2),'P(\Theta)')
hAx(2).YLim=[0 .1];

%Show Posterior
subplot(1,2,2)
plot(ThetaTest,PTheta_Data)
xlabel('\Theta')
ylabel('P(\Theta|Data)')

%% Markov-Chain Monte-Carlo

% web('https://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo')
% web('https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm')

N=10000;  %number of steps in  chain

%Intitialize chain vector
Chain=zeros(N,1);

%Setup initial value
Chain(1)=10;

%Setup probability function to explore
Mu=0;
Sigma=1;
p_x=@(x)normpdf(x,Mu,Sigma);

%Setup test jump size
TestSigma=.1;

%run chain

for nn=1:N-1
    %Make test jump
    CurrentX=Chain(nn);
    TestX=CurrentX+TestSigma*randn();
    
    %Metropolis-Hastings rule
    if rand()< (p_x(TestX)/p_x(CurrentX))
        Chain(nn+1)=TestX;
    else
        Chain(nn+1)=CurrentX;
    end
end

figure
hist(Chain,20)
xlabel('x')
ylabel('p(x)')

%% Visualize Chain

Xleft=-3;
Xright=11;

Xvec=(Xleft:.1:Xright);
Yvec=1:N;
FigH=figure;

for nn=N/2:N
    figure(FigH)
    %plot chain
    subplot(2,1,1)
    plot(Chain(1:nn),Yvec(1:nn))
    hold on
    plot(Chain(nn),Yvec(nn),'ro','Linewidth',3)
    axis([Xleft Xright 0 N])
    hold off
    xlabel('x')
    ylabel('Chain Index')
    
    %plot histogram
    subplot(2,1,2)
    [Counts]=hist(Chain(1:nn),Xvec);
    bar(Xvec,Counts)
    axis([Xleft Xright 0 max(Counts)*1.1])
    xlabel('x')
    ylabel('p(x)')
    
    pause(.001)
    
    
end





